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Abstract. The immersed- boundary method can be used to simulate flows 
around complex geometries within a Cartesian grid. This method has been 
used quite extensively in low Reynolds-numb er flows, and is now being 
applied to turbulent flows more frequently. The technique will be discussed, 
and three applications of the method will be presented, with increasing 
complexity, to illustrate the potential and limitations of the method, and 
some of the directions for future work. 


1. Introduction 

The increase in computer speed achieved over the last few’ years has made 
computational fluid dynamics increasingly useful and widespread as a tool 
to analyze and design flow’ configuration. Complex geometries, however, 
present an obstacle even to present-day computers, since the use of body- 
fitted meshes (structured or unstructured) significantly increases the cost 
of a calculation in terms of both computational speed and memory require- 
ments. 

An alternative method that may be cost-efficient in many situations is 
the fimmersed-boundary” method. This technique is based on the intro- 
duction of body forces distributed throughout the flow that mimic the effect 
that a solid body would have on the fluid. This approach allows the use 
of codes in Cartesian coordinates, which present significant advantages, in 
terms of speed, accuracy and flexibility, over codes that employ body fitted 
grids. 

This idea has been widely used in haemo-dynamics and bio-fluids en- 
gineering: two- and three-dimensional calculations of the flow’ in the heart 
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were reported by Peskin [13, 14] ajid McQueen and Peskin [8, 9]. In these 
calculations the motion of the boundary was determined by the fluid it- 
self, so that the boundary had to be modeled as a set of elements linked 
by springs. In cases in which the boundary motion is known a priori , the 
problem can be significantly simplified. 

Goldstein et a l. [4] proposed a feedback forcing mechanism that forces 
the fluid velocity u, to approach the velocity of the solid boundary. V), on 
the boundary itself. Consider the incompressible Navier-Stokes equations: 

(1) 

( 2 ) 

Goldstein et al. [4] assigned a force field 
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( 3 ) 

where a/ and j3j are two negative constants, and x Si j are the coordinates 
of the solid surface. The net effect of this force is to tend to annihilate 
the velocity difference u z — Vj. The flow, in fact, responds to the forcing 
as a damped oscillator (see [4] and [3] for an in-depth discussion of this 
issue); the frequency of the oscillator is oc \af\ 1 / 2 , whereas its damping 
coefficient is cc ,3// jal 1 / 2 . This implies that, in order to enforce the no-slip 
condition effectively, aj and 0f must have large magnitudes (larger than 
the highest frequency in the flow), which make the equations stiff. If the 
forcing is advanced in time explicitly, in fact, the maximum allowable CFL 
is o(10 -3 ); even implicit methods only allow the use of CFL numbers of 0.1 
or less. This constraint makes this method impractical for the calculation 
of unsteady (in particular, turbulent) flows. 

Recently, Mohd-Yusof [11] proposed the “direct forcing method/’ which 
assigns a force field given by 
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(4) 


(where the dependence on x S j and t has been omitted). This forcing im- 
poses directly the desired velocity on the immersed boundary, and has the 
advantage (over the feedback forcing method) that it does not require sig- 
nificant reductions in the allowable time-step. It was extensively tested in a 
staggered finite-difference code by Fadlun et al [3] , who derived an interpo- 
lation scheme to be used when the boundary does not coincide with a grid 
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Computed velocity ; 

outside the body 

Interpolated velocity 
....at the nearest outside pointL 



Figure. 1. Interpolation method used to apply the forcing. 


line. Verzicco et al [15] applied this method to the large-eddy simulation 
of high Reynolds number turbulent flows by calculating the flow inside an 
IC engine. 

In the present paper additional applications of this method will be pre- 
sented, and the potential and limitations of the technique, as well as issues 
that require further study, will be discussed. After a brief review of the 
governing equations and of the numerical scheme used, three test cases will 
be shown: a low Reynolds-number flow over a cylinder in the presence of a 
moving surface (Wannier [16] flow), the flow over a circular cylinder at low 
Reynolds number, and the bypass transition on a flat plate caused by the 
interaction between the boundary layer and the wake of a circular cylinder. 

2. Problem formulation 

Governing the flow are the incompressible Navier-Stokes equations (1-2). 
The flow solver is a standard 2nd-order accurate method on a staggered 
mesh [1]. The fractional time-step method [2, 6] is used and a 2nd-order 
accurate Adams-Bashforth method is employed for the time advancement. 
A non-reflecting boundary condition [12] is used at the outflow, and periodic 
boundary conditions in the spanwise direction. The inflow and free-stream 
conditions depended on the case studied. 

The direct forcing (4) was used. Since the immersed body does not 
follow a grid line the interpolation method proposed by Fadlun et al. [3] is 
used. The forcing is imposed not at the surface itself, but at the first point 
outside it (see Fig. 1) and the solid body velocity in (4) is replaced with 
the velocity Vi obtained by a linear interpolation between the computed 
fluid velocity outside the body, iq, and the desired body velocity Vj. This 
method has several advantages: first, it has been shown to be fully second- 
order accurate in time [3]; therefore, it is consistent with the second- order 
differencing scheme used by the solver; secondly, since it assumes that the 
velocity profile is linear near the body, it implies homogeneous Neumann 
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Figure 2. Wannier flow test case (Wannier 1950). (a) Computational domain and com- 
puted streamlines; (b) Li and L 2 norms of the error, e, for u and v velocity components. 
N is the total number of grid points. 

boundary conditions for the pressure (see the Appendix in [3: for a full 
discussion of this issue). This last feature is very important in the framework 
of fractional time-step methods, since it implies that the corrector step does 
not result in a modification of the body velocity imposed through the forcing 
in the Helmholtz step. On the other hand, the assumption that the velocity 
profile is linear over the first two layers of cells outside the immersed body 
requires the use of a very fine mesh in the vicinity of the body. 

3. Results and discussion 

3.1. WANNIER FLOW 

A straightforward way to verify the accuracy of the proposed methodology 
is to compute a flow containing a curved immersed boundary for which 
an analytical solution exists. The case considered here is the Stokes flow 
around a cylinder in the vicinity of a moving wall (see Fig. 2). An analytical 
solution for this case was derived by Wannier [16]. The streamlines for this 
flow are shown in Fig. 2a. Three computations on gradually finer uniform 
grids (32x32, 64x64, and 128 x 128) were conducted. The L\ and Lo norms 
of the error (the difference between the computed and analytical solution) 
are shown in Fig. 2b as a function of the total number of points N. The 
error decreases with a —2 slope indicating that the proposed methodology 
is second order accurate. 

3.2. FLOW OVER A CIRCULAR CYLINDER 

The next test case examined is the flow over a circular cylinder at Re& = 
UooD/u — 300 (where D is the cylinder diameter and Uqq the free- stream 
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Figure 3. Flow over a circular cylinder, Rep = 300. Velocity statistics, (a) U, (b) W, 
(c) (uV), (d) (u’w'). Reference data from Ref. [10]. 


velocity). Two calculations will be compared: a 2D one that used 400x200 
points in the xz— plane, and a 3D one, with the same mesh in the xz— plane, 
and 48 points in y. The computational domain was 60 x 27 t x 30. and the 
cylinder center was at x c = 10, z c = 15 (ail lengths are made dimensionless 
by £>, all velocities by Uqo). The grid was stretched both in the x — and 
z— directions; the last 1/6 of the domain (which required only 10 grid points 
in x) formed a sponge region, used to minimize the upstream propagation of 
disturbances due to the convective out Sow conditions. A uniform velocity 
profile was imposed at the inlet, and slip-wall conditions were applied at 


z = 0 and z = 30. 

The velocity statistics are shown in Fig. 3. Here and in the following the 
angle brackets denote averaging in time and in the spanwise direction. The 
3D calculation is in very good agreement with the reference data by Mittal 
and Balachandar [10]. At this Reynolds number, three-dimensionality is 
observed in the w T ake, which is evidenced in the visualization in Fig. 4, 
which shows iso- surfaces of the second invariant of the velocity- gradient 
tensor, 


1 dui duj 

2 dxj dx i 


2 \S\jSij HjjOvj) . 


( 5 ) 
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Figure 4 ■ Flow over a circular cylinder, Reo = 300. Isosurfaces of Q = 0.6. 


(where is the antisymmetric part of the velocity gradient tensor). The 
condition Q > 0 identifies effectively the regions of coherent vorticity [5] . 

Figure 4 shows the formation of an instability on the initially 2D rollers, 
and the presence of quasi-streamwise rib vortices joining the rollers. The 
magnitude of the spanwise Reynolds stresses { uV ). in this calculation, how- 
ever, remained significantly smaller than the other two normal components, 
which allowed the 2D calculation to give reasonable results. 

The effects of the grid resolution near the obstacle are shown in Fig. 5. 
The cell Reynolds number is defined as 

(6) 

where u and w are the instantaneous velocities, and Ax and Az the grid 
spacings. If the mesh is insufficiently fine ( Re c > 30), some oscillations 
can be observed that initiate along lines at ±45° on the cylinder (they 
are especially visible in the w contours, Fig. 5b). Refining the grid, thus 
reducing Re c reduces the size of this oscillation (Figs. 7a and b). Two- 
dimensional interpolation schemes that use both the points indicated by 
the diamonds and those indicated by squares in Fig. 6 to determine F; 
have also been found (Verzicco, private communication, 2001) to reduce 
the amplitude of these oscillations. 




Figure 6. One- dimensional vs. two-dimensional interpolation stencils. 

3.3. WAKE / B OUND ARY-LAYER INTERACTION 

Wakes interact with laminar boundary layers in many applications of en- 
gineering interest, for example on the leading edge of multi-component 




Figure. 8. Sketch of the wake/boundary-iayer configuration. 


airfoils, or inside turbo- machinery. The interaction of the turbulent eddies 
present in the wake with the boundary layer itself then becomes a primary 
driver of the transition process in the boundary layer itself, and may lead 
to transition to turbulence at fairly low Reynolds numbers. 

The configuration examined in the present study is sketched in Fig. 3. 
A circular cylinder, with its axis normal to the stream is placed above a 
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flat plate. The cylinder center is at x c = 10, z c — 3.2, immediately above 
the leading-edge of the plate, which was also at x c = 10. As in the previous 
case, distances are normalized by D, velocities by U oq. The computational 
domain was 60 x 27r x 20. As for the cylinder calculation, the grid was 
stretched both in the x — and z— directions and a sponge region was used. 
The Reynolds number based on cylinder diameter was 385. The configu- 
ration corresponds to Case 1 in the experimental paper by Kyriakides et 
al. [7], who observed significant velocity fluctuations in the boundary layer, 
starting from a point approximately six diameters downstream of the cylin- 
der. These fluctuations are generated by the large-scale convective motion 
of the vortices, and do not die down after the wake has weakened, but de- 
velop into a turbulent boundary layer despite the fact that the Reynolds 
number is very low. 

The distribution of the streamwise Reynolds stress, (u , u f ), as a function 
of x is shown in Fig. 9. A sudden increase of {v!v!) can be observed to begin 
at x — 8, indicating the beginning of transition. This result is consistent 
with the observations of Kyriakides et al. [7], who defined the onset of 
transition as “the :r— location where the velocity signal at the same height 
above the plate loses its sinusoidal character”, and found that transition 
occurs at x = 7.4. 

The velocity profiles, shown in Fig. 10a at several locations, initially 
resemble a Blasius profile merging into a wake near the cylinder. As the 
wake widens and interacts with the boundary layer, a logarithmic layer 
begins to establish itself, indicative of transition towards turbulent flow. 
This transitional behavior is also observed in the trace of the Reynolds 
stress tensor, q 2 (equal to twice the turbulent kinetic energy), which in 
the latter sections establishes a turbulent-like distribution, with a peak of 
magnitude 7 — 8^ at ~ 10 — 12. It should be noted that this quasi- 
turbulent state is achieved at very low Reynolds number: the boundarv- 



Figure 9. Wake/ boundary- layer interaction, {u'u') distribution at z = 0.13. 
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Figure 10. Wake/boundary-layer interaction. Turbulent statistics. Top: mean velocity 
profile; bottom: q = (u^ it'-). 


layer thickness 5 (defined as the distance above the plate at which the first 
maximum of the velocity profile occurs) is approximately 50-70 wail units. 

A visualization of the flow is shown in Fig. 11. The structure of the 
cylinder wake is similar to that highlighted in Fig. 4, with strong spanwise 
rollers that exhibit 3D instabilities and eventually break up, and smaller 
quasi-streamwise vortices in the braid region. The contours of streamwise 
velocity fluctuation u! on a plane parallel to the wall show significant; levels 
of fluctuations, especially near kinks in the rollers belonging to the lower 
row. Quasi-streamwise streaks are formed around x = 15, whose spacing is 
approximately 100 wall units. 

4. Conclusions 

The immersed- boundary technique has been presented and discussed. Il- 
lustrative results from three simulations indicate the potential of this tech- 
nique, which allows the calculation of flows around complex geometries 
without requiring a body-fitted grid. 
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If appropriate interpolation methods are used when the body does not 
coincide with a grid line, the method is second-order accurate. However, 
some care must be taken in the discretization of the 3ow field, especially in 
the vicinity of the body. We have observed the development of numerical os- 
cillations when the cell Reynolds- number exceeded values of approximately 
30. These oscillations did not grow or change position in time, and their 
effect on the flow field downstream of the obstacle was limited in the cases 
studied. This was observed, for instance, in calculations of the flow around 
the cylinder at Ren — 3500 were carried out in which the grid could not be 
sufficiently refined to satisfy the cell- Reynolds- number requirement. How- 
ever, it is not known whether these oscillations might give rise to instabili- 
ties in other geometries, or at higher Reynolds numbers. The development 
of more accurate, multi-dimensional interpolation schemes might be bene- 
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ficial in this respect. The use of multi-block methods, or embedded grids, 
could also alleviate this problem. 

If these numerical schemes can be overcome, the immericd midary 
method appears to be a useful tool for the simulation of flows in complex 
geometries at moderate or high Reynolds numbers. This is confirmed by 
the increasing number of studies using this method that are appearing in 
the literature. 
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